---
title: "Supplemental Information"
output: rmarkdown::null_document
---


```{r}

# Load Packages
library(pacman)

p_load(tidyverse, readxl, texreg, marginaleffects, ggpubr, stringr, scales, 
       Hmisc, broom, knitr, psych, xtable, kableExtra, ggrepel, ggforce, grid,
       ggsci)

# Student Sample -- SI Only
load("survey1_clean.rds")

# Gen Pop Sample -- Presented in Main Text 
load("survey2_clean.rds")



# Theme for Plotting 
my_theme <-  theme(
    text = element_text(size=14, color = "black"),
    strip.text = element_text(face = "bold"),
    axis.text.y = element_text(color = "black"),
    axis.text.x = element_text(color = "black"),
    axis.title = element_text(color = "black"),
    panel.background = element_blank(), 
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),  
    panel.border = element_rect(color = "black", fill = NA, linewidth = 1),  
    strip.background = element_rect(color = "black", fill = "white", linewidth = 1),
  )


```


# Descriptives 

## Make Table 2

```{r descriptive-by-condition}
survey2 <- survey2 |> 
  mutate(female = case_match(resp_gender,
                             "Female" ~ 1,
                             "Male" ~ 0),
         secondary = if_else(edu == "Post-Secondary", 1, 0),
         post_secondary = if_else(edu == "Secondary", 1, 0))


stats_s2 <- survey2 |> 
  drop_na(tx_scandal, tx_gender, edu, sel, 
          hostile, benevolent, left_right) |>
  group_by(tx_scandal) |> 
      summarise(
    obs = n(),
    across(c(hostile, benevolent, left_right, female,
             secondary, post_secondary), ~ round(mean(.x, na.rm = TRUE), 2))
  ) 


stats_s2 |> 
  pivot_longer(cols = obs:post_secondary,
               names_to = "index_type",
               values_to = "value") |> 
  pivot_wider(names_from = tx_scandal,
              values_from = value) |> 
  xtable(digits = 2) |> 
  print(include.rownames = FALSE, 
        file = "./tables/descriptives_s2.tex")
```


# Make Figure 1

```{r}

# Men take Sexual Initiative
init_stats <- survey2 |> 
  summarise(
    mean = round(mean(initiate_sex, na.rm = TRUE), 2),
    median = round(median(initiate_sex, na.rm = TRUE), 2),
    sd = round(sd(initiate_sex, na.rm = TRUE), 2)
  )

# Build plot
ggplot(survey2, aes(x = initiate_sex)) +
  geom_histogram(bins = 15, fill = "lightblue", color = "white") +
  labs(
    title = "“Normally, men are the ones who take the \ninitiative in sexual relationships.”",
    subtitle = "0 (Strongly Disagree) to 7 (Strongly Agree)",
    x = "Level of Agreement",
    y = "Number of Respondents",
    caption = paste0(
      "Mean = ", init_stats$mean, 
      " | Median = ", init_stats$median, 
      " | SD = ", init_stats$sd
    )
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold")
  )

ggsave("figures/si_plot_initiative.pdf", height = 4, width = 6)


```

# Conditional Effect of Hostile Sexism

## Estimate Models

```{r est-conditional-effects}

# Vector of DVs
dvs <- c("trust", "stay_involved", "vote")  

# Empty lists to store results for each scandal type
sex_results <- list()
corr_results <- list()
lie_results <- list()

# Loop through each DV and fit the model
for (dv in dvs) {
  
  # Create formula with each DV
  formula <- as.formula(paste(dv, "~ tx_gender * (benevolent + hostile + factor(edu) + sel + left_right)"))
  
  # Run model on each DV and store results separately
  sex_results[[dv]] <- lm(formula, data = subset(survey2, tx_scandal == "sexual-misconduct"))
  corr_results[[dv]] <- lm(formula, data = subset(survey2, tx_scandal == "corruption"))
  lie_results[[dv]] <- lm(formula, data = subset(survey2, tx_scandal == "lie"))

  }


```


## Estimate Slopes 

```{r est-slopes}
# Now estimate the average slopes using the results from above

# Dataframes for the results
sex_slopes <- tibble()
corr_slopes <- tibble()
lie_slopes <- tibble()

for (i in seq_along(dvs)) { # Looping over the DVs
  
  dv <- dvs[i]  # Get the dependent variable
  
  # Benevolent Sexism
  hostile_slopes_sex <- slopes(sex_results[[dv]], variables = "tx_gender",
                              newdata = datagrid(hostile = 0:7)) |> 
                      mutate(dv = dv, model = "sexual_misconduct") 
  
  
  hostile_slopes_corr <- slopes(corr_results[[dv]], variables = "tx_gender",
                               newdata = datagrid(hostile = 0:7)) |> 
                        mutate(dv = dv, model = "corruption")

  
  hostile_slopes_lie <- slopes(lie_results[[dv]], variables = "tx_gender",
                              newdata = datagrid(hostile = 0:7)) |> 
                        mutate(dv = dv, model = "lie")
  
  # Bind results together
  sex_slopes <- bind_rows(sex_slopes, hostile_slopes_sex)
  corr_slopes <- bind_rows(corr_slopes, hostile_slopes_corr)
  lie_slopes <- bind_rows(lie_slopes, hostile_slopes_lie)
  
}

# Bind Together Estimated Slopes
all_slopes <- bind_rows(sex_slopes, corr_slopes, lie_slopes)

# Label for plotting
all_slopes <- all_slopes |> 
  mutate(
    dv = case_match(dv, 
                    "trust" ~ "Trust Politician",
                    "stay_involved" ~ "Should Stay Involved", 
                    "vote" ~ "Vote for Politician"), 
    
    dv = factor(dv, levels = c("Trust Politician",
                               "Should Stay Involved", 
                               "Vote for Politician"))
    )



```


```{r plot-cond-effects}


# Function will Grab the Appropriate Scandal and Plot
slopes_plot <- function(scandal, title) {
  
  # Select the Scandal
  p_dat <-  all_slopes |> 
    filter(model == scandal) 
  
  # Plot 
  ggplot(p_dat, aes(x = hostile, y = estimate)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "lightblue", alpha = 0.3) +
  geom_line(color = "blue", size = 1) +
  labs(
      title = title,
      x = NULL,
      y = NULL,
    ) +
  facet_grid(~dv) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  theme_pubr() +
  theme(plot.title = element_text(hjust = .5))
    
}

# Combine Plots with ggarrange
combined <- ggarrange(
  slopes_plot(scandal = "sexual_misconduct", "Sexual Misconduct Scandal"),
  slopes_plot(scandal = "corruption", "Financial Corruption Scandal"),
  slopes_plot(scandal = "lie", "Dishonesty Scandal"),
  ncol = 1
)

# Add Labels 
annotate_figure(combined, bottom = text_grob("Hostile Sexism", 
                                             size = 14, face = "bold"),
                  left = text_grob("Male - Female", rot = 90, face ="bold", 
                                   size = 14))

# Save as Conditional Effects
ggsave("figures/conditional_effects_hostile.pdf", 
       height = 6,
       width = 8)


```


# Add Random State Fixed Effect to the Analysis

```{r add-res}


# Make Reference Category Nuevo Leon 
survey2 <- survey2 |> 
  mutate(state = relevel(state, ref = "Nuevo León"))

# Vector of DVs
dvs <- c("trust", "stay_involved", "vote")

# Empty lists to store results
main_models <- list() # list of main effects
interactive_models <- list()  # list of interactive effects

# Loop through each DV and fit the models
for (i in 1:length(dvs)) {
  dv <- dvs[i]

  # Main effects
  main_eff <- as.formula(paste(dv, "~ tx_scandal + tx_gender + factor(edu) + sel + hostile + benevolent + left_right + state"))
  main_mod <- lm(main_eff, data = survey2)

  # Interactive effects
  inter_eff <- as.formula(paste(dv, "~ tx_scandal * tx_gender + factor(edu) + sel + hostile + benevolent + left_right + state"))
  inter_mod <- lm(inter_eff, data = survey2)
  

  # Store results
  main_models[[dv]] <- main_mod
  interactive_models[[dv]] <- inter_mod

  
}


```



## Print to Latex

```{r print-tables-res}


# Main Effects
texreg(main_models,
      file = "tables/main_eff_w_state.tex",
      dcolumn = T,
      caption = "Main Effects including State Fixed Effects",
      caption.above = T,
      label = "si_main_eff_w_state",
      custom.model.names = c("Trust", "Stay Involved", "Vote"),
      custom.coef.map = list("(Intercept)" = "Intercept",
                              "tx_gendermale" = "Male Politician",
                              "tx_scandallie" = "Lie",
                              "tx_scandalsexual-misconduct" = "Sex",
                              "hostile" = "Hostile Sexism",
                              "benevolent" = "Benevolent Sexism",
                              "factor(edu)Primary/None" =  "Primary Edu.",
                              "factor(edu)Secondary" = "Secondary Edu.",
                              "sel" = "Socio-economic Level",
                              "left_right" = "LR Ideology"
                              ))



# Interactive Effects 
texreg(interactive_models,
       file = "tables/inter_eff_w_state.tex",
       dcolumn = T,
       custom.model.names = c("Trust", "Stay Involved", "Vote"),
       caption = "Interactive Effects including State Fixed Effects.",
       caption.above = T,
       label = "si_inter_eff_w_state",
       custom.coef.map = list("(Intercept)" = "Intercept",
                              "tx_gendermale" = "Male Politician",
                              "tx_scandallie" = "Lie",
                              "tx_scandallie:tx_gendermale" = "$\\times$ Male Politician", 
                              "tx_scandalsexual-misconduct" = "Sex Scandal",
                              "tx_scandalsexual-misconduct:tx_gendermale" = "$\\times$ Male Politican",
                              "hostile" = "Hostile Sexism",
                              "benevolent" = "Benevolent Sexism",
                              "factor(edu)Primary/None" =  "Primary Edu.",
                              "factor(edu)Secondary" = "Secondary Edu.",
                              "sel" = "Socio-economic Level",
                              "left_right" = "LR Ideology"))

```

## Plot the State Fixed Effects

```{r estimate-and-plot-fes}

# Empty DF for Results
state_fes <- data_frame()

# Loop Over DVS, Bind Together Estimates
for(i in 1:length(dvs)) {
  
  # State fixed effects comared to baseline 
  state_fe <- avg_comparisons(main_models[[i]], variables = list(state = "reference"))
  
  # Label the DV
  state_fe$dv <- dvs[i]
  
  # Bind Everything Together
  state_fes <- rbind(state_fes, state_fe)
  
}

# Clean State Labels for Plotting
plot_state_fes <- state_fes |> 
  mutate(State = str_remove(contrast, "- Nuevo León"),
         DV = str_replace(str_to_title(dv), "_", " ")) 

# Overall Mean of the Estimate, All Three DVs  
overall_mean_dv <- plot_state_fes |> 
  group_by(State) |> 
  summarise(overall_mean = mean(estimate)) 

# Join in Overall Mean
plot_state_fes <- plot_state_fes |> 
  left_join(overall_mean_dv)

# Plot the Fixed Effects
ggplot(plot_state_fes, aes(x = fct_reorder(State, overall_mean), 
                           y = estimate, 
                           color = DV,
                           ymin = conf.low,
                           ymax = conf.high)) +
  geom_pointrange() +
  #Facetiing by DV
  facet_row(~DV) +
  geom_hline(yintercept = 0, linetype = "dashed",
             color = "red", linewidth = .75) +
  coord_flip() + 
  labs(y = "Estimate", x = "State") +
  my_theme +
  theme(legend.position = "none")

ggsave(file = "figures/state_fes.pdf", height = 4.25, width = 6.5)
  

```



